using namespace std;

#include<iostream> 
#include<fstream>
#include<math.h>
#include<boost/property_tree/ptree.hpp>
#include<boost/property_tree/xml_parser.hpp>
#include<boost/foreach.hpp>

#include"flysim.h"
#ifdef USE_GPU
#include<CL/cl.hpp>
#include"opencl/clsafe.h"
#endif
#include"opencl/cl_constants.h"
#include"opencl/cl_interface.h"

// The following statements reserve address spaces for each of the 
// statically defined variables in cFlySim, the ancestor of all classes
// in the simulation.  They are set from a configuration file, which 
// also documents their usage.

static gsl_rng * rg = gsl_rng_alloc(gsl_rng_mt19937);
  
bool cFlySim::is_initialized = false;
// user specific variables
bool cFlySim::enable_gpu = false;
long cFlySim::ranseed = 0;
string cFlySim::outputfile;
string cFlySim::analysis;
int cFlySim::gciTotalNumberOfFlies;
int cFlySim::gciTotalNumberOfMales = 0;
int cFlySim::gciNumberOfPetriDishes = 0;
int cFlySim::debugSerialMPI = 0;
float cFlySim::dThreshold = 0;
float cFlySim::dMaleAveAve = 0;
float cFlySim::dMaleAveSD = 0;
float cFlySim::dMaleVarAve = 0;
float cFlySim::dMaleVarSD = 0;
float cFlySim::dFemaleAveAve = 0;
float cFlySim::dFemaleAveSD = 0;
float cFlySim::dFemaleVarAve = 0;
float cFlySim::dFemaleVarSD = 0;
float cFlySim::dMFCovarAve = 0;
float cFlySim::dMFCovarSD = 0;
float cFlySim::dMFCovCorrection = 0;

   //new user values
float cFlySim::dMpMAve = 0;
float cFlySim::dMpMSD = 0;
float cFlySim::dFpMAve = 0;
float cFlySim::dFpMSD = 0;
float cFlySim::dMpFAve = 0;
float cFlySim::dMpFSD = 0;
float cFlySim::dFpFAve = 0;
float cFlySim::dFpFSD = 0;
float cFlySim::dAvBoysWeight = 0;
float cFlySim::dAvGirlsWeight = 0;
float cFlySim::dVarBoysWeight = 0;
float cFlySim::dVarGirlsWeight = 0;
float cFlySim::dCovBoysGirlsWeight = 0;
float cFlySim::dMalesPerMaleWeight = 0;
float cFlySim::dFemalesPerMaleWeight = 0;
float cFlySim::dMalesPerFemaleWeight = 0;
float cFlySim::dFemalesPerFemaleWeight = 0;





int cFlySim::ciNoStats = 0; //the number of stats you want to store
float cFlySim::dSamplingInterval = 0;
int cFlySim::iNSamples = 0;
int cFlySim::iSamplingInterval = 0;
int cFlySim::iTimeToStartSampling = 0;
// ABC parameters
bool cFlySim::bDoABC = 0;
int cFlySim::gciNoOfIterationsInABCAnalysis = 0;
int cFlySim::gciBurnInPeriod = 0;
int cFlySim::gciABCSamplingPeriod = 0;
int cFlySim::gciNoOfSimulationsPerProposal = 0;
float cFlySim::dABCTolerance = 0;
float cFlySim::dABCNewTolerance = 0;
float cFlySim::dNewSumsSquaresValue = 0;
float cFlySim::dSquaredDistance = 0;
bool cFlySim::bUseBortotScheme = 0;
bool cFlySim::bToleranceChanged = 0;
float cFlySim::cdProbChangeTolerance = 0;
float cFlySim::cdProbToleranceMovesDown = 0;
float cFlySim::gdOldBortotTolerance; // book-keeps the value of the threshold when we are changing i = 0t
float cFlySim::gdMaxTolerance = 0;
float cFlySim::gdMinTolerance = 0;
float cFlySim::gdMaxToleranceChange = 0;
float cFlySim::gdInitialToleranceWhenUsingBortot = 0;
float cFlySim::gdPriorLambdaForTolerance = 0;
int cFlySim::gciNumberOfFlyParameters = 0;
bool cFlySim::bOutputSummaryStats = 0;
#ifdef USE_MPI
int cFlySim::mpi_rank = 0, cFlySim::mpi_ranks = 1;
#endif

const float cFlySim::SQRT_2 = sqrt(2);
//int cFlySim::iTotalTimeOfExperiment = 0;
gsl_rng * cFlySim::rg;
ofstream cFlySim::ofile;
cOpenCLInterface * cFlySim::clInterface;
//float * cFlySim::dPreComputedPowers;


float Random()
{
  //return gsl_rng_uniform_pos(rg);
  return 0;
}

cFlySim::cFlySim(){
  if (!is_initialized){
    throw "All flysim objects must be initialized before use\n";
  }
}

cFlySim::~cFlySim(){
  if (clInterface!=NULL) delete clInterface;
}


#ifdef USE_MPI
int cFlySim::get_rank(){
  return mpi_rank;
}
int cFlySim::get_ranks(){
  return mpi_ranks;
}

#endif

void cFlySim::initialize(const char * xmlfile){
        #ifdef USE_MPI
        if (MPI::COMM_WORLD.Get_rank()==0)
        #endif
  cerr<<"Initializing configuration from XML file\n";
  ifstream ifs(xmlfile);
  if (!ifs.is_open()){
    cerr<<"Could not open "<<xmlfile<<" configuration file\n";
    throw "Initialization failed\n";
  }
  ifs.close();
  clInterface = NULL;
    boost::property_tree::ptree pt;
    boost::property_tree::read_xml(xmlfile,pt);
    //load_from_xmlobject
    // load user values
    analysis = pt.get<string>("user.analysis");
 //   enable_gpu = pt.get<bool>("user.enable_gpu");
    ranseed = pt.get<long>("user.ranseed");
    outputfile = pt.get<string>("user.outputfile").data();
    gciTotalNumberOfFlies = pt.get<int>("user.gciTotalNumberOfFlies");
    gciTotalNumberOfMales = pt.get<int>("user.gciTotalNumberOfMales");
    gciNumberOfPetriDishes = pt.get<int>("user.gciNumberOfPetriDishes");
    debugSerialMPI = pt.get<int>("user.debugSerialMPI");
    dThreshold = pt.get<float>("user.dThreshold");
    dMaleAveAve = pt.get<float>("user.dMaleAveAve");
    dMaleAveSD = pt.get<float>("user.dMaleAveSD");
    dMaleVarAve = pt.get<float>("user.dMaleVarAve");
    dMaleVarSD = pt.get<float>("user.dMaleVarSD");
    dFemaleAveAve = pt.get<float>("user.dFemaleAveAve");
    dFemaleAveSD = pt.get<float>("user.dFemaleAveSD");
    dFemaleVarAve = pt.get<float>("user.dFemaleVarAve");
    dFemaleVarSD = pt.get<float>("user.dFemaleVarSD");
    dMFCovarAve = pt.get<float>("user.dMFCovarAve");
    dMFCovarSD = pt.get<float>("user.dMFCovarSD");
    dMFCovCorrection = pt.get<float>("user.dMFCovCorrection");
    //new user values
    dMpMAve = pt.get<float>("user.dMpMAve");
    dMpMSD = pt.get<float>("user.dMpMSD");
    dFpMAve = pt.get<float>("user.dFpMAve");
    dFpMSD = pt.get<float>("user.dFpMSD");
    dMpFAve = pt.get<float>("user.dMpFAve");
    dMpFSD = pt.get<float>("user.dMpFSD");
    dFpFAve = pt.get<float>("user.dFpFAve");
    dFpFSD = pt.get<float>("user.dFpFSD");
    dAvBoysWeight = pt.get<float>("user.dAvBoysWeight");
    dAvGirlsWeight = pt.get<float>("user.dAvGirlsWeight");
    dVarBoysWeight = pt.get<float>("user.dVarBoysWeight");
    dVarGirlsWeight = pt.get<float>("user.dVarGirlsWeight");
    dCovBoysGirlsWeight = pt.get<float>("user.dCovBoysGirlsWeight");
    dMalesPerMaleWeight = pt.get<float>("user.dMalesPerMaleWeight");
    dFemalesPerMaleWeight = pt.get<float>("user.dFemalesPerMaleWeight");
    dMalesPerFemaleWeight = pt.get<float>("user.dMalesPerFemaleWeight");
    dFemalesPerFemaleWeight = pt.get<float>("user.dFemalesPerFemaleWeight");

    ciNoStats = pt.get<int>("user.ciNoStats");
    dSamplingInterval = pt.get<float>("user.dSamplingInterval");
    iNSamples = pt.get<int>("user.iNSamples");
    iSamplingInterval = pt.get<int>("user.iSamplingInterval");
    iTimeToStartSampling = pt.get<int>("user.iTimeToStartSampling");
    // load ABC values
    bDoABC = pt.get<bool>("abc.bDoABC");
    gciNoOfIterationsInABCAnalysis = pt.get<int>("abc.gciNoOfIterationsInABCAnalysis");
    gciBurnInPeriod = pt.get<int>("abc.gciBurnInPeriod");
    gciABCSamplingPeriod = pt.get<int>("abc.gciABCSamplingPeriod");
    gciNoOfSimulationsPerProposal = pt.get<int>("abc.gciNoOfSimulationsPerProposal");
    dABCTolerance = pt.get<float>("abc.dABCTolerance");
    dABCNewTolerance = pt.get<float>("abc.dABCNewTolerance");
    dNewSumsSquaresValue = pt.get<float>("abc.dNewSumsSquaresValue");
    dSquaredDistance = pt.get<float>("abc.dSquaredDistance");
    bUseBortotScheme = pt.get<bool>("abc.bUseBortotScheme");
    bToleranceChanged = pt.get<bool>("abc.bToleranceChanged");
    cdProbChangeTolerance = pt.get<float>("abc.cdProbChangeTolerance");
    cdProbToleranceMovesDown = pt.get<float>("abc.cdProbToleranceMovesDown");
    gdMaxTolerance = pt.get<float>("abc.gdMaxTolerance");
    gdMinTolerance = pt.get<float>("abc.gdMinTolerance");
    gdMaxToleranceChange = pt.get<float>("abc.gdMaxToleranceChange");
    gdInitialToleranceWhenUsingBortot = pt.get<float>("abc.gdInitialToleranceWhenUsingBortot");
    gdPriorLambdaForTolerance = pt.get<float>("abc.gdPriorLambdaForTolerance");
    gciNumberOfFlyParameters = pt.get<int>("abc.gciNumberOfFlyParameters");
    bOutputSummaryStats = pt.get<bool>("abc.bOutputSummaryStats");
//	ciNumberOfGenerationsForSMCAnalysis=pt.get<bool>("smc.iNumberOfGenerationsForSMCAnalysis"); 
	//SMC::ciNoOfAcceptancesPerGeneration=pt.get<bool>("smc.ciNoOfAcceptancesPerGeneration");

#ifdef USE_GPU 
    if (enable_gpu){
      int platform_id = pt.get<int>("opencl.platform_id");
      int device_id = pt.get<int>("opencl.device_id");
      //cerr<<"Configuring with platform ID "<<platform_id<<" and device ID "<<device_id<<endl;
      string source_path = pt.get<string>("opencl.source_path");
      vector<string> source_vec;
      BOOST_FOREACH(boost::property_tree::ptree::value_type sources_object, pt.get_child("opencl.sources")){
        //boost::property_tree::ptree sources_tree = (boost::property_tree::ptree)sources_object.second;
        string source = sources_object.second.data();
        source_vec.push_back(source);
        //cerr<<"Added kernel source "<<source<<endl;
      }
      clInterface = new cOpenCLInterface(platform_id,device_id,source_path.data(),source_vec,ranseed);
      clInterface->test_rng();
    }
#endif
    // end load

  is_initialized = true;
}

/*float cFlySim::NormalRV(float mean,float sd){
  return gsl_ran_gaussian(rg, sd)+mean;
}*/

int main(int argc, char * argv[])
{
#ifdef USE_MPI
MPI::Init();
 // int mpi_rank = cFlySim::get_rank(), mpi_ranks = cFlySim::get_ranks(); // not sure why gary wants it this way..
int mpi_rank=MPI::COMM_WORLD.Get_rank();
int mpi_ranks=MPI::COMM_WORLD.Get_size();

#endif
  int arg=0;
  if (argc<2)
  {
    cerr<<"Usage: <XML configuration file>\n";
    return 1;
  }
  const char * configfile = argv[++arg];
  float dStartTime,dFinishTime;
  dStartTime=clock();  
  cExperiment * MyExpt = NULL;
  SMC * pSMC = NULL;
  try{
    cFlySim::initialize(configfile);
    if (cFlySim::analysis.compare("mcmc")==0){
      MyExpt = new cExperiment();
#ifdef USE_MPI
      if (mpi_rank==mpi_ranks-1)
#endif
      cerr<<"Proceeding with MCMC\n";
      MyExpt->run();
    }else if (cFlySim::analysis.compare("smc")==0 
//#ifdef USE_MPI         //FIXME 
//&& mpi_rank==mpi_ranks-1 
//#endif 
) {

#ifdef USE_MPI
        if (mpi_rank==0)
#endif
      cerr<<"Initializing SMC\n";
      pSMC = new SMC();
      pSMC->InitializeArrays();

#ifdef USE_MPI
	if (mpi_rank==0)
#endif
//      pSMC->PrintParameterValues(0);
      cerr<<"Proceeding with SMC\n";
      pSMC->RunSMCAnalysis();
      pSMC->PrintParameterValues(1);
    }else
#ifdef USE_MPI
if (mpi_rank==mpi_ranks-1)
#endif
{
      cerr<<"Invalid analysis mode of "<<cFlySim::analysis<<endl;
      throw "Can't initialize analysis";
    }
  }catch (const char * & e){
    cerr<<"Exception caught with message: "<<e<<endl;
  }catch (exception & e){
    cerr<<"General exception of "<<e.what();
  }

#ifdef USE_MPI
if (mpi_rank==mpi_ranks-1) {
#endif
  dFinishTime=clock();
  cerr<<"Run time: "<<(dFinishTime-dStartTime)/CLOCKS_PER_SEC<<"(seconds)"<<endl;
  if (MyExpt!=NULL) delete MyExpt;
  if (pSMC!=NULL) delete pSMC;
#ifdef USE_MPI
}
MPI::Finalize();
#endif
  return 0;

}


